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ABSTRACT 

We present a Bayesian approach to calculating the coefficients that convert the outputs of 
ALMA 183 GHz water- vapour radiometers into estimates of path fluctuations which can then 
be used to correct the observed interferometric visibilities. The key features of the approach 
are a simple, thin-layer, three-parameter model of the atmosphere; using the absolute measure- 
ments from the radiometers to constrain the model; priors to incorporate physical constraints 
and ancillary information; and a Markov Chain Monte Carlo characterisation of the posterior 
distribution including full distributions for the phase correction coefficients. The outcomes of 
the procedure are therefore estimates of the coefficients and their confidence intervals. We 
illustrate the technique with simulations showing some degeneracies that can arise and the 
importance of priors in tackling them. We then apply the technique to an hour-long test obser- 
vation at the Sub-Millimetre Array and find that the technique is stable and that, in this case, 
its performance is close to optimal. The modelling is described in detail in the appendices and 
all of the implementation source code is made publicly available under the GPL. 



1 INTRODUCTION 

The performance of millimetre and sub-millimetre wave interfer- 
ometers is often limited by the fluctuation of the properties of the 
Earth's troposphere along the lines of sight of each of the elements 
of the interferometer. If not corrected, these fluctuations lead to a 
fluctuating delay to each element and subsequent loss of correla- 
tion (and therefore sensitivity) and a limit on the maximum usable 
baseline length. Some recent simulations of the effect of these fluc- 
tuations on ALMA science have been presented by us (Nikolic et al. 
2008) and other authors (e.g., Asaki et al. 2005). 

ALMA plans to correct for these fluctuations by a combination 
of two techniques: 

(i) Fast-switching, that is interleaving science observations with 
observations of near-by phase calibrators that allow antenna phase 
errors to be solved for 

(ii) Water- vapour radiometry, that is observing the emission of 
atmospheric water vapour along the line of sight of each element of 
the array, and inferring and correcting from these observations the 
fluctuating path to each element. 

The current plan for ALMA is that fast-switching will be used 
with a cycle period of between around 10 and 200 seconds while 
fluctuations on timescales from 1 second up to the fast-switching 
timescale will be corrected by the water vapour radiometry tech- 
nique. The actual fast- switching time that will be used will depend 
on the weather conditions and the scientific requirements of each 
project as well as the achieved accuracy of phase correction with 
the water- vapour radiometers. Furthermore we expect to be able to 
use the WVRs to correct for the phase transfer errors, that is the 



errors that arises due to the different lines of sight to the phase cal- 
ibration and science sources. 

One of the key requirements for radiometric phase correction 
is a prescription for converting the observed sky brightnesses in the 
neighbourhood of the 183.3 GHz water vapour line, as measured 
by the WVRs, into a path delay that can be used to correct the 
observed visibilities. Any such prescription is complicated by the 
significant pressure-broadening and the high cross section of this 
line (it can be close to saturated even in the dry conditions of the 
ALMA site) which means that the optimum conversion strategy is 
quite a sensitive function of the prevalent atmospheric conditions. 

In this paper we assume that over the short timescales of inter- 
est (i.e., less than approximately 200 seconds) the differential path 
fluctuation between two telescopes can be predicted by a constant 
times the differential fluctuation of the sky brightness between the 
same telescopes. The constant of proportionality is the phase cor- 
rection "coefficient" and we give in this paper an initial prescription 
for calculating these coefficients (Section 2). We analyse this pre- 
scription with simulations (Section 3) and then apply it to a test 
observation at the Sub-Millimetre Array (Section 4). 



2 METHOD 

Our aim is to find estimates of, and confidence intervals on, the 
coefficients to be used for phase correction. We denote the coeffi- 
cients as dL/dTs i where L is the excess path to the telescope and 
Tq i are the sky brightnesses as measured by the four channels of 
the WVRs. We use this notation in general although this differential 
only makes mathematical sense when there is a known and invert- 
ible function which connect the cause of fluctuation in L with 7g 
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In this paper we assume that the fluctuation in path are caused by 
fluctuations of the total water vapour column only so that the fluc- 
tuation water vapour column is also the sole cause of fluctuations 
in T B j. 

The approach we take in this paper is to construct a physical 
model of the atmosphere with a number of unknown parameters 
and use some observables and the basic physical considerations to 
constrain the possible values the parameters. Once the distribution 
of possible values of model parameters are known, we can use the 
same model to compute the distribution of phase correction coef- 
ficients. The physical model employed in these initial studies is 
extremely simple: 

• We assume water vapour is the only cause of opacity and path 
fluctuation 

• We assume all the water is concentrated in a thin layer at a 
single pressure and temperature 

• We make the plane parallel approximation for computing ef- 
fects due to changes in elevation of the antenna 

This means that the model has only three unknown parameters, 
namely total water vapour column (c), temperature (T) and pres- 
sure (P). 

The reason for using such a simple model is that it allows very 
extensive computational analysis while containing the basic ingre- 
dients which we know must influence phase correction. As we gain 
experience with simulations and observations at the ALMA site, 
we intend to extend this model in directions which show an actual 
improvement in the inference of the phase correction coefficients. 

For the time being we will assume that the only true observ- 
ables that we have available to constrain the model are the four 
absolute sky brightness temperatures measured by the WVRs. In 
additional to the direct observables, we have some constraints on 
the possible values the model parameters can take; these are called 
priors. In principle the priors are a joint probability function of all 
of the model parameters, but in the present paper we simplify by 
taking them to be separable into a product of functions of one pa- 
rameter only, i.e., p(cTP) = p(c)p(T)p(P). 

The method of solution of this problem that we describe here 
is the standard Bayesian Markov Chain Monte Carlo approach. As 
usual, the inference problem is described by the Bayes equation 
(see, for example Jaynes 2003; Sivia & Skilling 2006): 



p(8\D) 



p(D\e)p(Q) 

PiP) ■ 



(1) 



where the symbols have following meaning: 



6 is the vector of model parameters (in this case {c, 7, P}) 

D is the observed data (in this case the sky brightness tempera- 
tures observed by the WVR) 

p(6) is the prior information (in this case constraints on model 
parameters as mentioned above) 

p(D\Q) is the likelihood, i.e., the probability of observing the data 
we have given some model parameters 6 

p(D) is the so-called Bayesian evidence, that is a measure of how 
good our model is in describing the data 
p(6\D) is the posterior distribution of model parameters 



2.1 Likelihood 

The computation of the likelihood of an observation given model 
parameters can be factored into three stages: 
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(a) Precipitable water vapour from 0.6 mm to 1.3 mm 
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(b) Temperature from 230 to 300 K 
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(c) Pressure from 400 to 750 mBar 

Figure 1. Variation of the shape of the 183 GHz water vapour line with 
changes in water vapour column, temperature and pressure 



(i) Calculation of the sky brightness temperature as a function 
of frequency 

(ii) Calculation of the temperatures recorded by the WVRs 
given a sky temperature. Here we model frequency response of the 
units and the coupling to the sky. 

(iii) Calculating the probability of observed data given the 
model temperatures obtained in the previous step 

The sky brightness is computed using a simplified model, in 
which we assume that the only relevant contributors to the atmo- 
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(b) Production radiometers 



Figure 2. Filter design centres and bandwidths of the prototype (top) and 
production (bottom) radiometers with the 1 83 GHz water vapour line also 
shown (in red). The heights of the rectangles representing the filters are 
inversely proportional to the square root of bandwidth, and are therefore an 
indication of the relative sensitivity of the filters. 



spheric opacity in this band are the water vapour line at 183 GHz 
and the water-vapour pseudo-continuum. We assume the water 
vapour line has a Gross line shape and parameters derived from 
the HITRAN database entry (Rothman et al. 2005) and suitable 
correction for pressure and temperature as detailed in Appendix A. 
The parametrisation that we use for the pseudo-continuum follows 
closely that in the program am by Paine (2004) and is also given 
in Appendix A. With the opacity calculation, the sky brightness 
can be calculated using simple radiative transfer. The resulting sky 
brightness temperature is shown for a range of conditions in Fig- 
ure 1 illustrating the change in the shape of the water vapour line 
with pressure, temperature and total water vapour column. 

The three model parameters are the total water vapour column 
(c), temperature (T) and pressure (P). The water vapour column 
parameter is taken to refer to the zenith column; for the cases when 
we investigate non-zenith measurements, we assume that this pa- 
rameter scales according to the plane-parallel approximation while 
the other parameters remain unchanged with a change in elevation 
of the telescope. 

With the sky brightness known, the temperature seen in each 
of the WVR filters can be calculated. In the present study we as- 
sume the WVRs are double- sideband (like the production units on 
ALMA will be) and that their filter set corresponds to the prototype 



filter set. The reason for using the prototype filter set definition is 
that later we will make an analysis of a sample test data set col- 
lected with the prototype radiometers at the SMA. The prototype 
filter set is shown in Figure 2, which also shows for contrast the 
filter set of the production ALMA WVRs. 

In this study we assume the filters are perfectly sharp and 
at their design centre frequencies and bandwidths. The average 
sky brightness across each filter sideband is calculated using the 
five point Gauss-Legendre Quadrature (e.g., Abramowitz & Stegun 
1964) which provides reasonable (but, since it is not an adaptive al- 
gorithm, a non-uniform) accuracy and is extremely efficient since 
the sky brightness needs to be calculated at only five frequencies 
per filter sideband. 

The second effect which needs to be taken into account at this 
stage is the non-perfect coupling of the radiometer in the sky. Our 
initial analyses of sky-dip measurements during this test campaign 
(to be published subsequently) suggest that the coupling was about 
0.91 and that the termination temperature of the parts of the beam 
that did not reach the sky was about 265 K. We therefore use these 
values in the analysis of the data from the SMA later (Section 4); 
however for the simulations shown in Section 3 for simplicity we 
assume perfect coupling. 

The results at stage of computation are the four temperatures, 
T& i, seen by the WVR as functions of 6 = {c,T,P}. Given these 
parameters 0, what is the probability of observing a set of temper- 
atures •? This is the probability that we denote by the likelihood 
function L and it is governed by the instrumental effects within the 
radiometer. 

For the present analysis, the dominant source of error is the un- 
certainty in the absolute calibration of the radiometers. The produc- 
tion radiometers specification require this error to be smaller than 
2 K at all times. If the underlying errors were normally-distributed 
(which may be a reasonable approximation) that would require an 
underlying a of less than 1 K. We do not however know at this time 
precisely what the final distribution will be and further it is likely 
that the calibration errors will be correlated between the filter out- 
puts. 

Nevertheless, for simplicity we presently take the errors to be 
normally and independently distributed so that the likelihood func- 
tion takes the well known form: 

logp(D|e) = logL = - £ ^ ^ j (2) 

where we take Cjy = 1 K. When we better understand the calibra- 
tion uncertainties in calibration of the WVR units, it will be im- 
portant to incorporate this information into the above equation. (By 
comparison the thermal noise in one second of integrating time will 
be below 0.1K.) 

2.2 Priors 

Besides the sky brightness observed by the WVRs, we have (and 
will have) some constraints on the model parameters that are the re- 
sults of physical considerations, or derived from independent past 
observations. For the purposes of this paper we will consider three 
different priors, shown in Table 1, with the aim of illustrating the 
effects they have on the final results. The priors are specific to the 
Mauna Kea site rather than ALMA site since they will also be used 
for analysis of testing data collected at the Mauna Kea. The conclu- 
sions derived from them are expected, however, to transfer directly 
to the ALMA site. 
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1 Basic p(c) = 

2 Reasonable p(c) = 

3 Pressure constraint p(c) = 



Omm < c < 5mm 
otherwise 

Omm < c < 5mm 
otherwise 

Omm < c < 5mm 
otherwise 



P(T) = 
P(T) = 
P(T) = 



200K<T< 320K 

otherwise 
260K<T<280K 

otherwise 
260K<T<280K 

otherwise 



P(P)- 
P(P)- 
P(P)- 



1 lOOmB < P < 650mB 

otherwise. 

1 530mB < P < 610mB 

otherwise. 

1 570mB < P < 590mB 
otherwise. 



Table 1. Priors used for the simulations and the analysis of SMA data 



The first prior we consider (number 1 in Table 1) is an ex- 
tremely relaxed prior which puts very loose constraints on the pres- 
sure and temperature of the water vapour layer: 

(i) Uniform probability that the pressure is between 100 and 
650mBar. Since the mean pressure at the peak of Mauna Kea is 
605 mBar this prior only assumes that the water vapour is not ex- 
tremely high in the atmosphere 

(ii) Uniform probability that the temperature is between 200 and 
320 K. Clearly this is a much wider range than typical tropospheric 
temperatures at altitudes where there is significant water vapour. 

(iii) Uniform probability that the zenith water vapour column is 
between and 5 mm. This uninformative constraint is used for all 
of the other priors too. 

Since this prior is less informational then the constraints we will 
have even without any ancillary measurements at the site, it is used 
to illustrate the degeneracies present if no priors at all are present. 

The second prior we consider (number 2 in Table 1) is de- 
signed to be representative of information on water vapour have 
might have with basic understanding of the site but without any 
sophisticated ancillary measurements. In this prior, we assume we 
know the temperature of the water vapour layer within 20 K and 
the pressure of the layer to within 80 mBar. This prior is used to 
analyse the test data from Mauna Kea. 

Finally, the third prior we consider (number 3 in Table 1) is 
used as an illustration of the improvement in accuracy that can be 
obtained with a tight prior on one of the parameters. In this case we 
still assume that we know the temperature to within 20 K but that 
we know the pressure of the water vapour layer to within 20 mBar 
instead of 80 mBar assumed in prior 2. Such a constraint on pres- 
sure of the water vapour layer might be derived from, for example, 
determination of its height using one of the techniques described 
below. 

2.3 Markov Chain Monte Carlo 

With the expressions for the likelihood and the priors that we have 
described above, we have the necessary information to compute 
the Bayes equation (Eq. 1). In general, this is computationally 
expensive (see for example, MacKay 2003) because of the large 
volume of parameter space that must be characterised in order to 
determine the maximum of p(D\0)p(6) and the numerical value 
p(D) = fd6p(D\6)p(d). 

The approach we take is standard Markov Chain Monte Carlo 
(MCMC) using the Metropolis et al. (1953) algorithm (for a tutorial 
review, see also Neal 1993). In this approach a chain of points in the 
parameter space is calculated such that the next point in the chain 
is found by proposing a new point by random displacement from 
the current point and calculating the relative likelihood of the two 
points. If the new point is more likely it is accepted onto the chain; 
if it is less likely, it accepted with probability determined by the 
ratio of the likelihoods of the current and new points. The proposal 
distributions we use in this paper are Gaussian with o c = 0.001 mm, 



Gt = 0.2 K and Gp = 0.5 mBar. We use the implementation of the 
MCMC algorithm in the BNMinl library (Nikolic 2009). 

By construction of the Metropolis algorithm, the density of 
points of the chain in a volume of parameter space is an estimate of 
the posterior probability p(6\D). This straight-forward approach 
does not however allow estimation of p(D) on its own, and so it 
is not possible to analyse the relative benefits of different models. 
There are however techniques available for estimation of p(D) and 
we intend to implement these in the future to allow proper model 
comparison. 

2.4 Marginalisation and calculating the dL/dTg 

We use the information contained in the Markov Chains in two 
ways: we marginalise and histogram the points to make estimates 
of the model parameters; and we, for each point in the chain, cal- 
culate the phase correction coefficients dL/dTg^. The marginalisa- 
tion of Markov Chains in directions parallel to parameter axes is 
trivially accomplished by simply ignoring those parameters. Sub- 
sequent computation of histograms is also easily done (we use the 
numpy . histogram routine for this). 

Computation of the coefficients dL/dTg^ is more involved as 
new physics must be introduced to compute the delay introduced 
by the water vapour layer. In the interest of simplicity, we split this 
calculation into two parts: calculation of the non-dispersive and the 
dispersive path delay. We compute the non-dispersive path delay 
using the Smith- Weintraub equation as described in Appendix B 1 
taking into account the temperature and the pressure of the water 
vapour. As we are interested here in the effects of water vapour 
only, we do not here consider the effect a change in temperature 
will have on the refractive index of the dry air. 

The dispersive delay calculation is more complicated as well 
as dependent on the observing frequency (unlike the rest of the dis- 
cussion presented in this memo). However, at the frequencies of 
relevance to the SMA test data presented later, the dispersive ef- 
fect is relatively small, i.e., around 5%. We therefore only calculate 
an adjustment using the ATM by Pardo et al. (2001) which is then 
used to scale the non-dispersive path. The details of this calculation 
are presented in Appendix B2 and the conclusion is that we scale 
up the non-dispersive phase coefficients by a factor of 1.05 to take 
into account the dispersion at 230 GHz. 

With the path delay calculated we calculate the coefficients 
dL/dTg i by making a small perturbation to the parameter c, i.e., 
the quantity of water vapour, and computing the differential as a 
finite difference. 



3 SIMULATION 

In this section we present analysis of a simulated single observation 
of the sky brightness with a WVR. The goals for this section are to 
illustrate the outputs of the technique we have described above, to 
show the effects and the importance of the priors, and to illustrate 
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the approximate accuracy with which it will be possible to predict 
the correction coefficients dL/dT^i from the inputs consisting of 
the sky brightness only. 

The atmosphere from which we simulate our data point has 
1 mm of water vapour toward the zenith at a temperature of 270 K 
and pressure of 580 mBar. We assume the observation is toward the 
zenith. As described above, we simulate the sky brightness mea- 
sured by the filter set of prototype WVRs. We find in this case that 
the simulated temperatures are: 194.8, 142.6, 90.7 and 47.5 K for 
the inner to outer channels respectively. These simulated sky tem- 
peratures are then used for subsequent inference as the observable 
temperatures (7^ •). 

The first inference we present is with the weak priors from row 
1 of Table 1. Recall that in this case we assume weak constraints on 
the model parameters, including allowing the pressure to be higher 
than the ground level pressure; hence, this case is an illustration 
of the inference when essentially no prior information is supplied. 
The posterior distribution of the model parameters for this case is 
shown in Figure 3. Like in the other figures of the model parameter 
posterior distribution, we present these results as the marginalised 
posterior distributions of each of c, T and P; and as the joint dis- 
tribution of c vs P, c vs T and c vs T with the remaining (third) 
parameter marginalised. 

The marginalised probabilities in the top row of Figure 3 show 
that in this case we can place relatively weak constraint on the 
amount of water vapour (to about 5%) and very poor constraints 
on both the pressure and the temperature of the water vapour. It can 
also be clearly noted that the posterior distributions are not well ap- 
proximated by a Gaussian distribution; for example, the posterior 
distribution p(c) shows a long tail toward higher water columns. 

The reason for the poor inference can be understood from the 
second row of Figure 3 which shows the joint distributions of each 
combination of the model parameters. Considering first the plot of 
the joint probability p(TP) in the right panel of the lower row, Fig- 
ure 3, we can see that the retrieval of the pressure and temperature 
are almost exactly degenerate. In other words if a certain combi- 
nation of pressure and temperature explain the observed data well, 
then a higher temperature and a proportionally higher pressure de- 
scribes the observation also sufficiently well. 

This degeneracy between pressure and temperature also af- 
fects the accuracy of the retrieval of the water vapour column as 
shown in the left and middle panels of the lower row of Figure 3. 
There we see that the extreme values of pressure and temperature 
that are permissible due to the degeneracy give rise to a tail of like- 
lihood to higher values of the retrieved water vapour column. 

The next result we describe combines the same simulated data 
point with priors in the second row of Table 1 . These are the priors 
that we might have without significant ancillary information, i.e., 
that we know the temperature of the water vapour layer to within 
20 K and its pressure to within 80 mBar. 

The posterior distribution from this inference is shown in Fig- 
ure 4. It can be seen from the top row of this figure that these re- 
sults are qualitatively different from the inference with very non- 
informative priors. In this case, the inference of the water vapour 
is well approximated by Gaussian with a full-width-half-maximum 
of about 0.012 mm and the entire distribution is within 0.02 mm of 
the model value. The inferences of the temperature and pressure are 
still poor however; in fact, it can be seen that their distributions fill 
almost entirely the space allowed by their priors, indicating that the 
priors in this case are providing important information. 

The lower row of Figure 4 again provides an explanation 
for the marginalised distributions of the model parameters. The 



pressure- temperature joint distribution again shows the degeneracy 
which explains the poor retrieval of each of those model param- 
eters individually. The water column-pressure and water column- 
temperature distributions still show spreads but with two important 
differences: 

(i) The priors mean that the range over which pressure and tem- 
perature can vary is much smaller, therefore leading to smaller er- 
rors in the retrieved water column 

(ii) The joint distributions are highly elongated along the ver- 
tical axes, which means a relatively large change in temperate or 
pressure is required to cause an error in the retrieved water column 

The condition (ii) is somewhat specific to the simulated point in 
parameter space, and will not be as true for a general combination 
of filter centres/widths and observed sky brightnesses. 

Since the inference shown in Figure 4 is constrained to a 
reasonable volume of the parameter space we can compute the 
dL/dTs i to find how well we can predict the phase correction co- 
efficients. The results are shown in Figure 5 as the marginalised 
distribution of each of the coefficients (upper part of the figure) 
and also the joint distribution of each pair of coefficients (lower 
part). It can be seen in the upper part of this figure that the pos- 
terior distributions are non-Gaussian, in fact almost square, and of 
width of about 10%. This large spread is probably dominated by the 
uncertainty in the retrieved temperature of the water vapour which 
influences its refractive index (see Equation B3). 

The joint distributions plotted in the lower part of Figure 5 
show that in this case the errors on inference of the coefficients 
are very highly correlated, with similar correlation for each pair of 
the channels. This means that it is unlikely that making use of all 
of the radiometer channels simultaneously would reduce the error 
in phase correction due to the uncertainty in the inference of the 
coefficients. Under different conditions and perhaps with additional 
observational data, this may be however be possible. 

The last result that we show in this section is an inference 
with a tight prior on the pressure (or, equivalently, the height) of 
the water vapour layer, i.e., 570mBar < P < 590mBar, but with 
the prior on temperature as before. The posterior distribution of the 
model parameters for this case is shown in Figure 6. 

We again find that the results of the inference are qualitatively 
different, primarily in that the posterior distribution of the temper- 
ature of the water vapour is now approximately Gaussian and with 
a full-width-half-maximum of about 8 K, significantly less than its 
prior range of 20 K. The posterior distribution for the pressure how- 
ever is, as expected, completely dominated by the prior and simply 
flat over its prior range. The implication is that a tighter a-priori 
range of the pressure allows a much better inference of the temper- 
ature of the water vapour layer. 

The posterior distributions of the dL/dT^ i coefficients for this 
posterior distribution are shown in Figure 7. The qualitatively dif- 
ferent inference of the temperature has an important impact on the 
quality of the inference of these coefficients too, as can be seen 
from the approximately Gaussian shape of coefficients for chan- 
nels 1 to 3 and by their significantly tighter ranges. For example, the 
distribution of coefficient for channel 2 for the inference with prior 
2, (Figure 5) is essentially flat with a range of about 0.008 mm/K. 
With a tight prior on the pressure however, the distribution of this 
coefficient is close to Gaussian with a full-width-half-maximum of 
0.003 mm/K, clearly substantially better. 
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Figure 4. Posterior distribution of model parameters derived from four absolute sky brightness temperatures only, i.e., like Figure 3, but now including 
reasonable flat priors on all three of the c, T and P model parameters. 



3.1 Discussion 

One of the main aims of this section was to illustrate the outputs 
of an analysis based on the methods described in Section 2, i.e., 
viewing the problem of estimating the phase correction coefficients 
as Bayesian inference problem. 

The main outcome of such an analysis are the posterior dis- 
tributions for the phase correction coefficients such as presented in 
Figures 5 and 7. The posterior distributions allow us both to pick 
a specific coefficient to use for each channel of the radiometer and 
give us confidence intervals for the accuracy of those coefficients. 
Obtaining such confidence intervals is important since the combi- 



nation of dynamic scheduling, wide range of ALMA configurations 
and a large number of projects will mean that each science obser- 
vation is likely to be in conditions which are just marginal for its 
requirements. Therefore, if phase correction doesn't work as well 
as expected, it is likely to seriously impact the science aims. 

We should note however that the confidence intervals calcu- 
lated using the model we have specified and therefore do not prop- 
erly capture the probability that simply our model is not very good. 
Within the Bayesian framework this can be done through the evi- 
dence, or p(D), value and we intend to implement this functionality 
in the near future. 
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(b) Join distributions of the phase correction coefficiens 

Figure 5. Posterior distributions of the coefficients dL/dTg^ used to convert brightness fluctuations into path fluctuations for the model parameters posterior 
shown in Figure 4. 



Beside the distributions of the phase correction coefficients, 
the outcome of the analysis is also the full joint distribution of the 
model parameters, which in this case are the water vapour column, 
and its pressure and temperature. It is the availability of this full 
posterior distribution that makes reliable estimates of the phase co- 



efficients possible. This becomes particularly important when more 
parameters are added to the problem, as will no doubt be necessary 
in our case. With more parameters it becomes less and less feasi- 
ble to pick a single representative point in the parameter space to 
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Figure 6. Model parameter posterior distribution with tight prior on pressure. 
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Figure 7. The posterior distribution of (XL/dT^ i corresponding to Figure 6. 



calculate the phase correction coefficients at and making use of the 
full distribution becomes increasingly important. 

The model we have been using in this paper is fairly simplis- 
tic in that it assumes that the only observables we have are the four 
absolute sky brightness temperatures observed by the radiometers. 
What we find is that if we have no further constraints at all, we can 
estimate the water vapour column to an accuracy of about 5%. We 



can not however constrain the temperature of the water vapour at 
all because it is almost exactly degenerate with the pressure and 
therefore we can not make any estimate of the phase correction co- 
efficients (since they depend on the temperature, see Equation B3). 

This shortcoming can be improved on by adding even a fairly 
loose constraint on the pressure and temperature of the water, such 
as can be derived from historical distributions of water vapour in 
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the atmosphere and its temperature. Such loose constraints should 
be able to provide estimates of the phase correction coefficients in 
the 10% range. Even tighter a-priori constraints can provide further 
improvements as shown in Figures 6 and 7. 

With the model as simple as the one we have presented here 
and no data from the site of ALMA itself it would be somewhat 
premature to discuss in detail already the improvements in accuracy 
particular constraints can make. We however note that we will have 
a substantial amount of ancillary information such as: 

• Ground level air temperature, pressure and relative humidity 
and a number of points at the site 

• Inference of the vertical temperature profile of the atmosphere 
at the centre of the array from a commercial O2 line sounder 

• Library of vertical profiles of atmospheric parameters from 
past radiosonde launches 

• Meso- scale meteorological forecasts 

• Inference of atmospheric parameters from specialised tele- 
scope observations such as sky-dips and crossed beam observations 

Inferences from these measurements can be used as a-priori prob- 
abilities on model parameters, or, indeed, in some cases, as further 
observables which are analysed simultaneously with the sky bright- 
ness measurement. 

It should be noted however that further and better prior infor- 
mation on the temperature and pressure parameters will leave two 
other sources of uncertainty: errors due to the inaccuracies in the 
model that we use and limitations in estimating the water vapour 
column which arises due to calibration accuracy of the WVRs. The 
best way of tackling these uncertainties may turn out to be to use 
the observed correlation between changes in the excess path to the 
telescopes (8L) and the fluctuations of the temperature observed 
by the radiometers (STq i) as an additional observable that can con- 
strain the models. We will discuss this approach in the next memo 
in this series. 

Finally, we now consider ways in which the model presented 
above may be improved. Firstly, the model may fairly easily be 
re-parametrised in terms of the height of the water vapour layer, 
the temperature lapse rate of the atmosphere and the ground level 
pressures and temperatures instead of the current parametrisa- 
tion in terms of pressure and temperature of the layer itself. This 
would allow us to easily take into account the measured ground- 
level temperature and pressure and the further information that 
we will have on the temperature lapse rate (through historical 
radiosonde measurements and oxygen line profiling) and water 
vapour height (through historical radiosonde data and specialised 
telescope scans). The second improvement is to consider the effect 
of a thick layer of water vapour, perhaps with an exponential fall 
off in water vapour content as a function of height. 

Assessing the benefits of such improvements to the model will 
require computation of the evidence value and, ideally, test data 
from the ALMA site. 



4 APPLICATION TO TEST DATA COLLECTED AT SMA 

In this section we analyse an observation taken with the two pro- 
totype ALMA water vapour radiometers (Hills et al. 2001) at the 
Submillimetre Array (SMA) on Mauna Kea. 

The observation was made on 17 February 2006 and con- 
sisted of an approximately one hour long track of a bright quasar. 
The interferometric visibility between the two SMA antennas with 
the water- vapour radiometers was recorded together with the sky 
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Figure 9. Retrieval of the zenith water vapour column as a function of time 
from the data recorded by one of the radiometers during the February 17th 
observation. The priors used in the retrieval are the priors from row 2 of 
Table 1 . The retrieval was calculated once every 25 seconds for the hour- 
long observation. 



brightnesses seen by the radiometers. This observation was part 
of the testing campaign of the prototype radiometers at the SMA 
which will be described in detail in a separate paper. 

The effective LO frequency of the observation was 235 GHz 
and both the upper and the lower sideband were recorded. In the 
present analysis we use only the upper side band with the on-sky 
frequency of 240 GHz. The water- vapour radiometers were read out 
with an integration time of 1 second while the interferometer was 
read out at 2.5 seconds. In order to match the two data sets, the 
radiometer data ware re- sampled to 2.5 second time resolution, and 
a small adjustment to the time-stamps was made to correct for a 
known timing drift problem. The radiometer data were recorded in 
already calibrated format and required no further processing. The 
interferometric data were converted to a text format by M. Reid of 
the SMA. Subsequently the variation in the observed visibility was 
transformed to the path fluctuation between the two antennas. 

The total sky brightness temperatures observed by four chan- 
nels of two radiometers during this observation are shown in Fig- 
ure 8. It can be seen that observed brightness is increasing during 
the observation which is a consequence of the decreasing eleva- 
tion of the source and therefore increasing airmass as the observa- 
tion progressed. It can also be seen that the innermost two channels 
(channels 1 and 2, top row of Figure 8) were almost saturated, indi- 
cating significant water vapour along the line of sight. The overall 
levels of the blue and red curves in this plot allow us to make infer- 
ence about the atmospheric conditions at the time of the observation 
while the relative fluctuations between the two curves, once multi- 
plied by the phase correction coefficients, should correlate closely 
with the path fluctuation measured by the interferometric observa- 
tion. 

A retrieval of atmospheric parameters can be made from ev- 
ery one second integration of each of the two radiometers without 
significant loss in accuracy since, as mentioned before, the thermal 
noise in one second is much smaller than the expected calibration 
error of the radiometers. In practice however we expect retrievals 
will be made rather less often since it is likely that variation in at- 
mospheric conditions will be fairly slow at the ALMA site. In this 
study we will present retrievals at three different time -resolutions: 
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Figure 8. The sky brightness temperatures observed by the two prototype WVRs while at tests at the SMA, which was tracking a quasar at the time. Blue is 
one of the radiometers and red is the other. The four panels represent the four radiometer channels. 



(i) Calculation of the marginalised zenith water vapour at 300 
time points during the observation 

(ii) A detailed analysis of a single data point in the middle of the 
observation 

(iii) Calculation of the marginalised phase correction coeffi- 
cients at three points of observation 

We first present a sequence of 300 retrievals from sky bright- 
ness temperatures measured at intervals spaced by 25 seconds dur- 
ing the observation. As with the other retrievals in this section, 
these were made with the priors as shown in row 2 of Table 1. 
We plotted these retrievals in Figure 9 by marginalising them to 
obtain the posterior distribution of the zenith water vapour column 
and plotting these histograms in colour scale as a function of time. 
Therefore in this plot, time runs in the horizontal direction, the 
zenith water vapour column parameter along the vertical and the 
colour represents the posterior probability. 

The first point to note is that because we retrieve for the zenith 
water vapour column, we do not expect to see a large increase in the 
parameter c as the air-mass increases during the observation. The 
overall range of the water vapour measured in the retrievals is about 
5% indicating that at least approximately the plane-parallel approx- 
imation holds and the referencing to zenith column is a reasonable 
approach. 

Secondly it can be noted that although fluctuations in the wa- 
ter vapour column can be seen, adjacent retrievals (which are sep- 
arated by 25 seconds in time) generally show similar values in- 
dicating a high degree of stability in the retrieved posterior dis- 
tribution. The magnitude of the fluctuations seen in the retrieved 
water column is about 50 ^im of water vapour which corresponds 
approximately to a 300 \im of path fluctuation (using the rough con- 



version of 1 mm water « 6 mm path). As we will present later, e.g., 
Figure 13, this roughly corresponds to the fluctuations seen in the 
path by the interferometric measurements. Therefore, over the hour 
of the observation and a large change in airmass there is no obvious 
evidence of instability in the water vapour column retrieval. 

A more detailed analysis of the retrieval from one set of sky 
brightness measurements in the middle of the observation is pre- 
sented in Figure 10 in the same format as the previous plots of the 
posterior distributions of model parameters. As the line of sight 
water vapour during this observations is significantly higher than 
the simulations shown in Section 3, some qualitatively differently 
different results are seen. 

The first noticeable feature is that in this retrieval the temper- 
ature of the water vapour layer is in fact quite well constrained, to 
a range of around 4 Kelvin, in contrast the results at 1 mm water 
vapour line of sight column where the posterior distribution filled 
the entire prior range of 20 K. This is a direct consequence of the al- 
most complete saturation of the innermost channel of the radiome- 
ter which means that it is in effect measuring the temperature of the 
water vapour rather than line of sight column. 

The saturation of the inner channel, however, now also causes 
a degeneracy between the pressure and the water vapour column pa- 
rameter leading to a non-Gaussian posterior distribution for c. The 
reason for this is that at lower pressure, the line opacity becomes 
highly peaked, but this can not be detected because of the satura- 
tion. The outcome therefore is that at this point in parameter space 
the water vapour column is somewhat less well constrained and the 
temperature is better constrained. It should again be noted that the 
inference errors we present are derived from the model itself and 
therefore do not account for inaccuracies of the model. 
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We next consider the posterior distribution of the phase correc- 
tion coefficients for this detailed retrieval, as shown in Figure 1 1 . 
Because of the high line of sight water vapour column, it can be 
seen that the coefficient for channel 1 is very high, i.e., around 1.1 
mm/K. This means that thermal noise of 0. 1 K within the radiome- 
ter would be sufficient to produce a high random path fluctuation 



of 110 \im. Furthermore, other effects which are not modelled are 
likely to cause large errors in the path fluctuations calculated from 
this channel. Therefore, in this case we do not expect this channel 
to provide useful data for the phase correction itself. 

It can also be seen in Figure 1 1 that although the posterior dis- 
tributions of the phase correction coefficients are reasonable cen- 
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trally peaked, they do all have a tail to higher values. This tail is a 
result of the non-Guassian inference of the total water vapour col- 
umn, as seen in the top-left panel of Figure 10. The errors on these 
coefficients are not dominated by the error on inferred water vapour 
temperature because, as discussed above, the saturated inner chan- 
nel in allows for this error to be minimised in this case. Finally it 
should be noted that the confidence interval of the inferred phase 
correction coefficients is about 10 to 15% in this case. 

The final part set of inferences we made was to repeat the pre- 
diction of the phase correction coefficients at the beginning and at 
the end observation so that we can examine how they change dur- 
ing the course of the observation. The results of this analysis are 
shown in Figure 12, with each retrieval on a separate row. 

It can be seen in that the values of all of the coefficients in- 
crease as the observation progresses. This is of course due to the in- 
crease in airmass. The greatest increase is seen in channel 1, which 
starts at a median value of about 0.45 and finishes at a median value 
of about 2.25. 

4.1 Phase correction 

We next examine how well we would have been able to do phase 
correction at the SMA during this observations if we had used the 
phase correction coefficients retrievals as shown previously. 

As mentioned previously the SMA interferometer made it pos- 
sible to infer the fluctuation of path between the two antennas of the 
array that the WVRs mounted on them. We can also difference the 
signals between two radiometers to calculate the differential fluc- 
tuation of sky brightness as seen from the two antennas. Those two 
signals are compared (for each of the four channels of the WVRs) 
in Figure 13. Of course we hope that when multiplied by the phase 
correction coefficient, the sky brightness fluctuations will look very 
similar to the path fluctuation. 

The comparison shown in Figure 1 3 can immediately confirm 
that channel 1 is unlikely to yield a useful phase correction, while 
channel 2 show some correlation with path, and channels 3 and 4 
show very high correlation. 

We next quantify the performance that the phase correction 
technique would have produced based on the phase correction co- 
efficients inferred in this section. For this we consider both the ob- 
servation as a whole and three sub-sets of the observation, namely 
the first, the middle and the last thirds. The reason for considering 
the subsets is to investigate if the change in the phase correction 
coefficients inferred during the observation (Figure 12) is reflected 
in improved correction performance. 

Before processing the data further we remove from both the 
interferometric path and the sky temperature measurements a three- 
minute running mean. This trend removal should approximately 
model the effect of the further offset calibration scheme that ALMA 
will use in practice. Although we do not present them here, the data 
without the running mean removal lead to generally the same con- 
clusions as below. 

We next asses the potential performance of the phase correc- 
tion. For these tests we consider phase correction using the fluctua- 
tions observed in only one channel at the time. In principle, using a 
combinations of channels rather than a single one should give bet- 
ter performance because of the reduced effect of thermal noise and 
possibly the averaging out of errors in the inferred phase correction 
coefficients. We do pursue this multiple-channel approach further 
in this paper but we expect it will be used in practice with ALMA. 

We calculated the corrected path fluctuation using the inferred 
phase correction coefficients and we also find the 'best-fit' phase 



coefficient, i.e., the coefficient which would have produced smallest 
corrected path fluctuation given the data that we have. The results 
are present in Table 2 with following columns: 

(i) The root-mean- square of the uncorrected path fluctuation as 
seen by the interferometer, <7 raw 

(ii) The approximate median inferences for the phase correction 
coefficients for each channel and for each observation subset 

as shown in Figure 12 

(iii) The root-mean- square path fluctuation after correction us- 
ing radiometric data from each channel, o Y i 

(iv) The phase correction coefficient that would have given the 
optimum reduction of path fluctuation, Opt( ) 

(v) The resultant root-mean- square path fluctuation that would 
have been obtained if the optimum correction coefficient were used, 

°opt,i 

As can be seen from the table, channel 1 can not be used in 
these conditions used for phase correction, as expected from the re- 
sults already presented. This is evident in phase correction with the 
'optimum', i.e., best-fit, phase correction coefficient which makes 
a negligible improvement to the path fluctuation, reducing it from 
157 micron to 156 ]im. If the phase correction were done with the 
inferred coefficient of 1.25 mm/K (this was calculated for the mid- 
dle of the observation) then the path fluctuation would have been 
drastically increased. This is due to the high effect of thermal noise 
in this channel when the line is essentially saturated and other ef- 
fects such as temperature variations of the water vapour. 

The results for the other three channels show that significant 
improvements in corrected versus uncorrected path fluctuations can 
be made based on the radiometer data. Furthermore, it can be seen 
path corrected with our inferred coefficients dL/dT^ i is close to the 
path corrected with the best-fit coefficient, i.e., is fairly close to the 
limit that be achieved with the present data. The difference between 
using the inferred and best-fit coefficients results in changes of as 
little as 1 \im rms in path for some of the data sections up to 7 ^im 
path rms in the worst combination of section and filter. 

Some trends can however be noticed in the data presented in 
the table. Firstly, excluding channel 1 which is saturated, the in- 
ferred phase correction coefficients are systematically underesti- 
mating the best-fit coefficients. The only exception is for channel 
2 in the last third of the observation, which may however be also 
affected by saturation. Secondly there is evidence that at lower el- 
evations the gap between correction using inferred and best-fit co- 
efficients widens, suggesting that refinements to the model may be 
necessary. 



4.2 Discussion 

The analysis presented in Table 2 shows that the Bayesian infer- 
ence approach that we described in Section 2 can produce very 
useful phase correction coefficients even with a very simple atmo- 
spheric model. We have shown that the inferred coefficients give a 
corrected path fluctuation to within about five percent of what the 
optimal coefficient would have given. 

There are however a number of interesting points that these 
data (and the simulations) raise. Firstly, the confidence interval for 
the inference from our model analysis (Figure 12) show a range of 
about 10% and furthermore the optimal coefficients found empir- 
ically are about 15% higher than the median inferred values. This 
needs to be compared against the specification for ALMA phase 
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Figure 12. Change in the posterior distribution of the four phase-correction coefficients during the one-hour long scan tracking a quasar. Top row is beg ginning 
of the scan (UT 17.1), middle row is the middle of the scan (UT 17.5) and bottom row is the end of the scan (UT 17.8). 
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(b) Part II 

Table 2. Results of analysis of SMA data from February 17th. The part I of the table shows results for channels 1 and 2; the second part for channels 3 and 4. 
The UT range column indicates the time range of the data analysed: the first row corresponds to the whole set and the other three rows to respective thirds of 
the data. The column cr raw is the RMS of phase fluctuations without correction. The column is the best estimate from retrieval of the correction coefficient 
for the ith channel and for that section in time (see also Figure). The column 0" r i is the RMS phase fluctuation after correction using the estimated coefficient. 
The column Opt( ) is the coefficient which would give the minimum RMS give the radiometer and interferometer data. The column 0" opt i is that optimum 
RMS. 



correction which is: 

^corrected < f 1+ T~~i 10^m + 0.02 X 8Lr aw . (3) 

If we interpret the first term on the right hand side as the budget for 
the thermal noise contribution to corrected phase fluctuations then 
the second term, which is proportional to the raw phase fluctuation, 
among other effects (e.g., Nikolic et al. 2007) must account for the 
imperfect phase correction due to errors on the inferred phase cor- 
rection coefficients. Therefore we need to be able to infer phase cor- 



rection coefficients to better to 2% accuracy, but the present analy- 
sis of the model and the data suggest that we will do significantly 
worse than that. 

There are a number of improvements that we are planning that 
will allow better inference of the phase correction coefficients. One 
of these is the use of the observed correlation between the change 
in path and sky temperature which is likely to allow a significant 
improvement of the inference accuracy and will be presented in 
one of the next papers in this series. Secondly, it should be noted 
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Figure 13. Comparison of the differences between each of the channels (parts a-d) on the two radiometers at SMA test with the fluctuation of path observed 
by the interferometer. 



that in this case one of the limiting factors on the inference is the 
absolute calibration of the radiometers. Although we do not think 
that the a-priori calibration of the radiometers when mounted on 
the ALMA telescopes will be better than what we assumed here, it 
will be probably the case that the major sources of calibration inac- 



curacy will be fixed in time. Therefore we may be able to improve 
the accuracy of these devices over time by empirical corrections to 
their calibration. 

We have not discussed in detail in this paper the dispersive 
path delay, in particular how it depends on the parameters other 



2009 ALMA Memo 587 



15 



than the water vapour column. This can be investigated by more 
fully featured atmospheric models along the lines of the simple cal- 
culations performed in Appendix B2. We have also not considered 
the 'dry' path fluctuations, that is path fluctuation due to changes 
in refractive index of dry air. It is not presently clear how signif- 
icant in practice these will be but data to be collected during the 
commissioning of ALMA should provide constraints on this. 



5 CONCLUSIONS 

We have described an approach for calculation of the phase correc- 
tion coefficients based on: 

• A very simple, essentially minimal, model of the atmosphere 
with three parameters: water vapour column, pressure and temper- 
ature 

• Observables which are the four absolute sky brightness tem- 
peratures measured by the water vapour radiometers 

• Priors on the model parameters based either on basic consid- 
erations or ancillary measurements 

• Bayesian inference 

The approach produces the full probability distributions of the 
phase correction coefficients making it easier to asses in a timely 
way the accuracy of phase correction that can be expected. 

The simulations we presented in Section 3 show that al- 
though there are inherent degeneracies in retrieval of atmospheric 
paramours from sky brightness only, with the appropriate use of 
prior information fairly good retrievals may be obtained. 

We have then applied this approach to one test observation 
taken during the testing of the ALMA prototype radiometers at the 
SMA. We found that the inferred water vapour column during this 
observation is stable and its fluctuations are not greater than ex- 
pected given the observed path fluctuations (Figure 9). In tests of 
applying the phase correction, we found that good phase correc- 
tions can be made by taking representative inferences of the phase 
correction coefficients and applying them to correct the observed 
path fluctuation: in fact the typical corrected path fluctuation is 
within about 5% of the best correction that could be achieved with 
these data, assuming only one channel is used at a time. In order for 
it to be possible to meet the very demanding ALMA specifications, 
it is clear to we will need to make use of further data and possibly 
more complicated models. The approach we are currently working 
on is to incorporate information from the occasional phase calibra- 
tion observations of the observed correlation between path and sky 
brightness fluctuations into the inference process. 

The source code for all of the algorithms presented in this 
paper is available for public download under the GNU Pub- 
lic License at http : //www . mrao . cam . ac . uk/~bn204/alma/ 
memo-infer.html. The observational data from the SMA is also 
available at the same internet address. 
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APPENDIX A: CALCULATION OF ATMOSPHERIC 
OPACITY 

Al The 183 GHz water vapour line 

The water line at 183 GHz is assumed to have the Gross line shape: 



K(v;S,v ,y) : 



V 2 ) 2 +4 7 2v2' 



(Al) 



where v is the frequency and the set are {S, Vo, y} are the parame- 
ters of the line. 

In the present model we are assuming that water vapour is a 
relatively small constituent of air and so that the volume mixing 
ratio for water is negligible. In this case, the line parameters can be 
computed as follow: 



v = F + PxA a] 

?ref ^ 



7=7air XPX 



W- 7 **) l-e~ hv / T Q(T Kf ) 
Hxe r/ ret- - lLU 



l- e -hv/T M Q(J) • 



(A2) 
(A3) 

(A4) 



Here T re f is the reference temperature of the catalogue, which in this 
case is 296 K; T and P are the temperature and pressure of the air 
containing the water vapour, and Q(T) is the partition function of 
water vapour which we take from a look-up table. The remaining 
parameters are taken from the HITRAN entry for the 183.3 GHz 
water vapour line and are shown in Table Al. 



Of the remaining two terms in Equation B 1 , the dominant term 
is the second and it can be transformed to be in terms of water 
vapour density (see Stirling et al. 2008): 

yR 



dn ? 



(B2) 



where Jt v = 18.02gmol _1 is the molecular weight of water vapour 
and p w is the density of the water vapour. Assuming a thin layer and 
by integrating along the line of sight we get: 

1741 K 



61 1 



yR 



dc 



-de. 



(B3) 



B2 Dispersive path delay 

The dispersive delay adjustment was calculated using the ATM code 
Pardo et al. (2001). The calculation was made using a simple pro- 
gram dispersive that takes a frequency and ground-level relative 
humidity and computes the dispersive and non-dispersive delay due 
to the water vapour and the dry atmosphere at the supplied fre- 
quency. 

The program dispersive, and the code for the library 
ATM, are available at http: //www.mrao . cam. ac .uk/~bn204/ 
alma/ at momo del .html, packaged as aatm. All of the code is 
made available under the GPL license. For this paper we used 
aatm-0 . 06, i.e., version 0.06. 



A2 Continuum 

The continuum cross section is : 

C(v;Cb,r ,m)=Cb^) xAfxv 2 , (A5) 

where N is number density of air, and the parameter set {Co, 7b, m} 
describe the continuum. For water, the values we assume are: 
Param Value Unit 

C 6.1 xl0~ 48 cm 5 /GHz 2 
T 300.0 K 

m 2.6 none 



APPENDIX B: EXCESS PATH LENGTH 



Bl Non-dispersive path delay 

The non-dispersive excess path length due to water is computed 
using the Smith-Weintraub equation which gives the following ex- 
pression for the excess refractive index (Stirling et al. 2008): 



1 = 10" 



A 



-P- 



'ft 



(Bl) 



where P$ is the partial pressure of the dry air and P w is the partial 
pressure of the water vapour and the constants have the following 
values: 

• a = 11.6 x lO^KPa- 1 

• p =64.8x 10- 2 KPa _1 

• y= 3.776 x 10 3 K 2 Pa" 1 . 

We are presently concerned with the excess path introduced with 
the water vapour only, so we can drop the first term of Equation B 1 . 
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Parameter 



Units 



Value 



Comment 



H 

F 
Aair 

Ex 

Tair 
/self 
"aii- 



Cm 2 x GHz 
GHz 
GHz x mbar -1 
K 

GHz x mbar -1 
GHz x mbar -1 
none 



2.333884 x 10 -21 

183.310107 
-5.0298 x 10 -5 

195.908398 
2.9114 x 10 -3 
1.6149 x 10 -2 
0.77 



Line intensity 
Frequency 
Air-induced freq. shift 
Lower state energy level 
Air broadening 
Self-broadening 
Temperature dep. of air broadening 



Table Al. The HITRAN parameters and values for the 183.3 GHz line (Rothman et al. 2005). 
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